#!/usr/bin/perl -w
use strict;
use warnings;
use lib "/net/cpp-group/Leo/bin";
use parse_bl2seq;
use parse_fasta;
use Getopt::Long;

#originally
#parse_codeml_positions.pl		Scott Beatson 		20.05.03
# modified by Leo Goodstadt
#
#A program to parse the output from multiple codeml mlc files to generate
#a fasta formatted file that can be used as a footer in an alignment program

my $usage = <<'USAGE';

USAGE:

parse_slr_results.pl 
				slr.results
				--significant_residues output.fa
				> results.summary
				[--verbose]
				[--help]

    Parses the output from Slr results and prints
    a fasta formatted file indicating +ve selected positions.

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $verbose = undef;
my $significant_residues;
GetOptions(
			'help'			   		=> \$help,
			'verbose'		   		=> \$verbose,
			'significant_residues=s'	=> \$significant_residues,
			);

die $usage if ($help);
die $usage unless ($significant_residues);
die $usage unless (defined($ARGV[0]));

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print positive residues as FASTA
#
sub print_fasta_line(\@$$$$)
{
	my ($sites, $name, $cutoff, $file, $correct) = @_;
	print $file ">$name selected residues $correct p < $cutoff\n";
	my ($cnt_sites) =0;
	for my $site_index(0.. @$sites - 1)
	{
#		print $file $cnt_sites++,"\t", ($site < $cutoff ? '*' : '-') , "\n";;
		print $file ($sites->[$site_index] < $cutoff ? '*' : '-');
	}
	print $file "\n";


}
sub print_significant_residues(\@\@\@\@)
{
	my ($neg_sites, $pos_sites, $neg_sites_corrected, $pos_sites_corrected) = @_;
	open (SITES, ">$significant_residues") or
		die "Error:\n\tCould not open the file for saving positively / negatively selected residues in FASTA format [$significant_residues]\n$!\n";

	for my $cutoff(0.05, 0.01, 0.001)
	{
		print_fasta_line(@$neg_sites, "Negatively", $cutoff, *SITES, "");
	}
	for my $cutoff(0.05, 0.01, 0.001)
	{
		print_fasta_line(@$pos_sites, "Positively", $cutoff, *SITES, "");
	}

	for my $cutoff(0.05, 0.01, 0.001)
	{
		print_fasta_line(@$neg_sites_corrected, "Negatively", $cutoff, *SITES, "corrected");
	}
	for my $cutoff(0.05, 0.01, 0.001)
	{
		print_fasta_line(@$pos_sites_corrected, "Positively", $cutoff, *SITES, "corrected");
	}


}
#
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



my $ignore_first_line = <>;
my $cnt_sites=0;
my @sites_pos;
my @sites_neg;
my @sites_pos_corrected;
my @sites_neg_corrected;
my (
	$cnt_negative_050,  
	$cnt_negative_010, 
	$cnt_negative_001, 
	$cnt_positive_050,  
	$cnt_positive_010, 
	$cnt_positive_001,
	$cnt_negative_050_corrected,  
	$cnt_negative_010_corrected, 
	$cnt_negative_001_corrected, 
	$cnt_positive_050_corrected,  
	$cnt_positive_010_corrected, 
	$cnt_positive_001_corrected,
	)  = (0,0,0,0,0,0,0,0,0,0,0,0);
while (<>)
{
	
	my (undef, $site_num, undef, undef, $omega, undef, $pval, $pval_corrected, undef) = split /\s+/;
	die unless $site_num == ++$cnt_sites;
	if ($omega  < 1.0)
	{
		if ($pval_corrected < 0.05)
		{
			$cnt_negative_050_corrected++;
			if ($pval_corrected < 0.01)
			{
				$cnt_negative_010_corrected++;
				if ($pval_corrected < 0.001)
				{
					$cnt_negative_001_corrected++;
				}
			}
		}
		push(@sites_neg_corrected, $pval_corrected);
		push(@sites_pos_corrected, 1.0);


		if ($pval < 0.05)
		{
			$cnt_negative_050++;
			if ($pval < 0.01)
			{
				$cnt_negative_010++;
				if ($pval < 0.001)
				{
					$cnt_negative_001++;
				}
			}
		}
		push(@sites_neg, $pval);
		push(@sites_pos, 1.0);
	}
	else
	{
		if ($pval_corrected < 0.05)
		{
			$cnt_positive_050_corrected++;
			if ($pval_corrected < 0.01)
			{
				$cnt_positive_010_corrected++;
				if ($pval_corrected < 0.001)
				{
					$cnt_positive_001_corrected++;
				}
			}
		}
		push(@sites_pos_corrected, $pval_corrected);
		push(@sites_neg_corrected, 1.0);

		if ($pval < 0.05)
		{
			$cnt_positive_050++;
			if ($pval < 0.01)
			{
				$cnt_positive_010++;
				if ($pval < 0.001)
				{
					$cnt_positive_001++;
				}
			}
		}
		push(@sites_pos, $pval);
		push(@sites_neg, 1.0);
	}
}



print "\t$cnt_negative_050\tsites negatively selected at p < 0.05\n";
print "\t$cnt_negative_010\tsites negatively selected at p < 0.01\n";
print "\t$cnt_negative_001\tsites negatively selected at p < 0.001\n";
print "\t$cnt_positive_050\tsites positively selected at p < 0.05\n";
print "\t$cnt_positive_010\tsites positively selected at p < 0.01\n";
print "\t$cnt_positive_001\tsites positively selected at p < 0.001\n";
print "\t$cnt_negative_050_corrected\tsites negatively selected at corrected p < 0.05\n";
print "\t$cnt_negative_010_corrected\tsites negatively selected at corrected p < 0.01\n";
print "\t$cnt_negative_001_corrected\tsites negatively selected at corrected p < 0.001\n";
print "\t$cnt_positive_050_corrected\tsites positively selected at corrected p < 0.05\n";
print "\t$cnt_positive_010_corrected\tsites positively selected at corrected p < 0.01\n";
print "\t$cnt_positive_001_corrected\tsites positively selected at corrected p < 0.001\n";

print_significant_residues(@sites_neg, @sites_pos, @sites_neg_corrected, @sites_pos_corrected);

